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Abstract 

We  present  a  (suboptimal)  filtering  algorithm  for  tracking  a  highly  maneuvering  target  in  a 
cluttered  environment  using  multiple  sensors  dealing  with  possibly  asynchronous  (time  delayed) 
measurements.  The  filtering  algorithm  is  developed  by  applying  the  basic  Interacting  Multiple 
Model  (IMM)  approach,  the  Probabilistic  Data  Association  (PDA)  technique,  and  asynchronous 
measurement  updating  for  state-augmented  system  estimation  for  the  target.  A  state  augmented 
approach  is  developed  to  estimate  the  time  delay  between  local  and  remote  sensors.  A  multi¬ 
sensor  probabilistic  data  association  filter  is  developed  for  parallel  sensor  processing  for  target 
tracking  under  clutter.  The  algorithm  is  illustrated  via  a  highly  maneuvering  target  tracking 
simulation  example  where  two  sensors,  a  radar  and  an  infrared  sensor,  are  used.  Compared 
with  an  existing  IMMPDA  filtering  algorithm  with  the  assumption  of  synchronous  (no  delay) 
measurements  sensor  processing,  the  proposed  algorithm  achieves  considerable  improvement 
(especially  in  the  case  of  larger  delays)  in  the  accuracy  of  track  estimation. 

Keywords:  Asynchronous  (Delayed)  Measurements;  Multisensor  Parallel  Updating;  Interact¬ 
ing  Multiple  Model  (IMM);  Probabilistic  Data  Association  (PDA). 

I  Introduction 

We  consider  the  problem  of  tracking  a  single  maneuvering  target  in  clutter.  This  class  of  problem 
has  received  considerable  attention  in  the  literature  [1,  2,  3,  4,  9].  In  target  tracking  systems 
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measurements  axe  typically  collected  in  “scans”  or  “frames”  and  then  transmitted  to  a  processing 
center  [5,  6].  Asynchronous  (delayed)  measurements  arise  in  a  multisensor  central  tracking  system 
due  to  communication  network  delays,  varying  preprocessing  times  at  the  sensor  platforms  and 
possibly  lack  of  sampling  time  synchronization  among  sensor  platforms.  One  of  the  asynchronous 
measurement  problems  is  that  of  out-of-sequence  measurements  (OOSM)  where  measurements  at 
various  sensors  may  arrive  out-of-sequence  (not  in  correct  time  order)  at  the  central  processor. 
OOSM  has  been  considered  using  interacting  multiple  model  (IMM)  [6,  7,  8].  In  this  paper  we  do 
not  consider  OOSM  but,  instead,  consider  “in-sequence”  measurements  with  a  fixed  but  unknown 
relative  time-delay  among  sensor  measurements.  Various  sensor  measurements  are  assumed  to  be  at 
the  same  rate  but  not  necessarily  time  synchronized.  All  measurements  over  one  sampling  interval 
(based  on  the  local  clock  of  the  central  processor)  are  collected  at  the  central  processor,  attributed 
to  one  time  instant  and  processed  simultaneously.  We  exploit  interacting  multiple  model  (IMM) 
and  probabilistic  data  association  (PDA)  techniques.  It  is  assumed  that  a  track  has  been  formed 
(initiated)  and  the  objective  of  this  work  is  to  investigate  fixed-but- unknown  relative  time-delay 
(measurement  timing  mismatch)  arising  in  a  multisensor  central  tracking  system. 

In  [6],  fixed- lag  smoothing  techniques  have  been  investigated  using  IMM  algorithm  combined 
with  PDA  filter  in  a  multiple  sensor  scenario  to  propose  a  combined  IMMMSPDAF  (interacting 
multiple  model  multiple  sensor  probabilistic  data  association  filter) .  We  exploit  the  basic  structure 
of  [1]  in  combination  with  a  state-augmented  approach  to  deal  with  the  fixed-but-unknown  relative 
time-delay.  In  [1]  and  [14]  it  is  assumed  that  the  sensors  are  collocated  and  (time)  synchronized 
with  the  sampling  rate.  In  contrast,  the  sensor  collocation  and  (time)  synchronization  are  no  longer 
assumed  in  this  paper.  Also,  unlike  [1,  9,  12]  which  have  used  sequential  updating  of  the  state 
estimates  with  measurements  (i.e.,  updating  of  the  state  estimates  sequentially  with  measurements 
from  different  sensors),  we  use  parallel  updating  of  the  state  estimates  with  measurements  (i.e., 
updating  of  the  state  estimates  with  all  measurements  at  the  same  time).  For  linear  systems,  the 
two  updating  methods  are  algebraically  equivalent  but  for  nonlinear  filtering,  the  parallel  updating 
can  yield  better  performance  in  spite  of  higher  computational  cost  [4].  Ref.  [14]  uses  parallel 
updating  but  has  some  errors:  during  data  association,  all  measurements  at  the  same  time  from 
different  sensors  are  assumed  to  be  either  from  clutter  or  from  the  target.  The  possibility  that 
a  measurement  from  sesnsor  1  may  be  from  target  while  the  measurement  from  sensor  2  may  be 
clutter-induced  (and  vice-versa)  is  implicitly  not  allowed  in  [14]  -  this  is  clearly  incorrect.  Ref.  [10] 
allows  for  such  distinctions  (hypotheses),  however,  it  is  limited  to  non-manevering  targets.  In  this 
paper,  we  also  extend  the  multisensor  approach  of  [10]  to  maneuvering  targets  (see  Step  4  in  Sec. 
IV). 

The  paper  is  organized  as  follows.  Section  II  presents  the  problem  formulation.  Section  III 
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describes  the  state-augmented  system  approach.  Section  IV  describes  the  proposed  IMMMSPDAF 
algorithm  for  asynchronous  measurements.  Simulation  results  using  the  proposed  algorithm  for  a 
realistic  problem  are  given  in  Section  V.  Finally,  Section  VI  presents  a  discussion  of  the  results 
and  some  conclusions. 

II  Problem  Formulation 

We  assume  that  the  target  dynamics  can  be  modeled  by  one  of  n  hypothesized  models.  The  model 
set  is  denoted  as  Mn  :=  {1, ...,  n}  and  there  axe  total  q  sensors.  The  event  that  model  m  is  in  effect 
during  the  sampling  period  (tk-i,tk]  is  denoted  by  M™.  For  the  mth  hypothesized  model  (mode), 
the  state  dynamics  and  measurements,  respectively,  are  modeled  as 

=  *SVi**-i  +  G™k- lC-i  (1) 

and 

4  =  hm'l(xk)  +  w™’1  for  l  =  :  local  model  at  the  sensor,  (2) 

where  xk  is  the  system  state  at  tk  and  of  dimension  nx,  zlk  is  the  (true)  measurement  vector  (i.e.,  due 
to  the  target)  at  sensor  l  at  tk  and  of  dimension  n2;,  Fk^k_1  and  G1kk_l  are  the  system  matrices  when 
model  m  is  in  effect  over  the  sampling  period  (tfc-i,  ijJ,  and  hm’1  is  the  nonlinear  transformation  of 
xk  to  zlk  (/  =  1,  ...,g)  for  model  m.  A  first-order  linearized  version  of  (2)  is  given  by 

4  =  H?'lxk  +  w™'1  for  l  =  (3) 

where  H™'1  is  the  Jacobian  matrix  of  hm’1  evaluated  at  some  value  of  the  estimate  of  state  xk 
(see  Sec.  III).  The  process  noise  and  the  measurement  noise  w™'1  are  mutually  uncorrelated 
zero-mean  white  Gaussian  processes  with  covariance  matrices  Qkl_1  and  R™’1,  respectively.  At  the 
initial  time  to,  the  initial  conditions  for  the  system  state  under  each  model  m  are  assumed  to  be 
Gaussian  random  variables  with  the  known  mean  E{x™}  and  the  known  covariance  P™.  The 
probability  of  model  m  at  to,  n™  —  P{M™},  is  also  known.  The  switching  from  model  Ml_x  to 
model  M™  is  governed  by  a  finite-state  stationary  Markov  chain  with  known  transition  probabilities 
Pim  =  P{ilC|MLi}.  Henceforth,  time  tk  will  be  denoted  by  k. 

Assume  that  there  is  a  fixed  but  unknown  relative  time  delay  dk  (modulo  resampling  interval) 
at  sample  time  tk  between  the  local  sensor  clock  and  the  central  processor  clock  at  sample  time 
tk.  [This  time  delay  could  be  due  to  unsynchronized  clocks  at  the  two  locations  or  due  to  inherent 
delay  due  to  congestion,  insufficient  bandwidth  etc.  in  the  communication  link  between  the  remote 
sensor  platform  and  the  central  processor.]  The  measurements  from  sensor  l  are  sent  to  the  central 
processor  where  all  measurements  collected  between  local  sampling  interval  (tfc_i,  tk]  are  attributed 
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to  time  tk.  The  state  dynamics  and  measurements  reported  from  the  remote  sensor  platform  at 
time  tkdl  (henceforth  will  be  denoted  by  kdi)  to  the  center  processor  at  time  tk  can  be  modeled  as 


xkdi  -  F&^k-iZk-i  +  Gkdl,k-i» k-\  (4) 

and 


4  =  hm,l(xkdl)  +  w™’1  :  model  of  the  sensor  at  the  central  processor  (5) 

where  tkdl  =  tk—dki  and  dki  is  the  time  difference  between  the  sampling  time  at  the  central  processor 
and  the  measurement  time  at  the  local  sensor  (assume  that  0  <  dki  <T,  where  T  is  sampling  time), 
Xkdl  is  the  system  state  at  tkdl  and  of  dimension  »*,  F an<^  ^ku,k-i  8X6  system  matrices 
when  model  m  is  in  effect  over  the  timing  interval  (tk-utkdl} 

III  State- Augmented  System 

Define  the  augmented  state  xk  from  xk  as 

4  =  [4>4>4-i>4-i]  (6) 

where  x'k  denotes  the  transpose  of  xk.  Assume  that  there  is  a  fixed  but  unknown  delay,  dki ,  between 
the  central  processor  and  the  remote  sensor  l  platform.  Using  the  above  definitions  (1,  6)  and  the 
measurement  delay,  dki,  the  augmented  state  equation  may  be  written  more  compactly  as 

^  =  (7) 

and 

dki  =  (8) 

where  is  a  small  processing  noise  assumed  to  be  Gaussian  noise  with  zero  mean  and  (very) 
small  but  nonzero  variance.  Note  that  the  process  noise  in  (7)  is  v™  (at  time  k  not  at  time  k  —  1). 
Above  equations  (7)  and  (8)  can  also  be  absorbed  into  another  augmented  state  xk  as 


Xk 


xk 

dki 


=  Fk!k-ixk~i + &k<k- i«r  where  = 


„dl 

Jk-1 


(9) 


x'k  =  [x'k,v'k,x'k_x,v'k_x ,dki],  and  Fk,k-i  and  Gk,k- 1  are  defined  in  Sec.  V  (see  (46)-(53)).  Using  the 
augmented  state  (9)  the  counterparts  to  (2)  and  (5),  respectively,  are 

4  =  hm’l(xk)  +  w™’1  =  hm'l([1, 0, 0, 0, 0]xfc)  +  w™'1  (10) 

and 
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for  both  measurements  from  local  sensor  and  from  remote  sensor,  respectively.  To  keep  the  notations 
and  details  to  a  bare  minimum,  we  will  consider  the  case  of  two  sensors  only  and  furthermore,  we 
will  assume  that  one  of  the  sensors  is  either  collocated  with  or  is  synchronized  with  the  central 
processor,  so  that  we  will  drop  the  subscript  l  from  dki-  For  more  than  two  sensors,  we  need  to 
augment  ik  with  additional  <4’s  (total  q  —  1):  in  essence,  these  delays  are  relative  to  one  of  the 
sensors  (reference  sensor). 

The  following  notations  and  definitions  are  used  regarding  the  measurements  at  sensor  l.  Note 
that,  in  general,  at  any  time  some  measurements  may  be  due  to  clutter  and  some  due  to  the  target, 
i.e.  there  can  be  more  than  a  single  measurement  at  time  k  at  sensor  /.  The  measurement  set  (not 
yet  validated)  generated  by  sensor  l  at  time  k  is  denoted  as 

4~{  4m,4m . 4<m,))  02) 


where  mi  is  the  number  of  measurements  generated  by  sensor  l  at  time  k.  Variable  z1^  (i  = 
l,...,mj)  is  the  ith  measurement  within  this  set.  The  validated  set  of  measurements  of  sensor  l  at 
time  k  will  be  denoted  by  Tjf,  containing  fhi  (<  m{)  measurement  vectors.  The  cumulative  set  of 


validated  measurements  from  sensor  l  up  to  time  k  is  denoted  as 

ZKI)  :={Y(,Y‘,...,Yl}.  (13) 

The  cumulative  set  of  validated  measurements  from  all  sensors  up  to  time  k  is  denoted  as 

Zk:={ZkW,Zk(-2\...,Zk («)}  (14) 

where  q  is  the  number  of  sensors. 

Our  goal  is  to  find  the  state  estimate 

%k\k  E{ik\Z  }  (15) 

and  the  associated  error  covariance  matrix 

Pk\k  -=  ®fc|fc][^/c  ®/c|fc]  \Z  }  (1®) 

where  x'k  denotes  the  transpose  of  a?*,.  * 


IV  IMMMSPDAF  Algorithm  for  Asynchronous  Measurements 

We  now  modify  the  IMM/(J)PDA  algorithms  of  [9]  and  [12]  to  apply  to  the  multi-sensor  asyn¬ 
chronous  measurements  system.  We  confine  our  attention  to  the  case  of  2  sensors;  however,  the 
algorithm  can  be  easily  adapted  to  the  case  of  arbitrary  q  sensors.  We  will  only  briefly  outline  the 
basic  steps  in  “one  cycle”  (i.e.,  processing  needed  to  update  for  a  new  set  of  measurements)  of  the 
IMMMSPDA  filter. 
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Assumed  available :  Given  the  state  estimate  :=  E{xk-\\M™_X,  Zk~1},  the  associated 

covariance  and  the  conditional  mode  probability  :=  P[M^l1|Zfc-1]  at  time  fc  —  1  for 

each  mode  m  6  Mn. 

Step  1.  Interaction  -  mixing  of  the  estimate  from  the  previous  time  (Vm  6  Mn)  : 
predicted  mode  probability: 

:=  mnz*-1]  =  Ewu-  (w) 

i 

mixing  probability: 


n"m  :=  P{MU\M£',Zt-']  = 


mixed  estimate: 


*EVi  := 


covariance  of  the  mixed  estimate: 

j£Vl  :=  -  itut-ilfe-l  -  iEVll'IMT. z‘_1> 

t 

Step  2.  Predicted  state  and  measurements  for  sensors  1  and  2  (Vm  €  Mn) 
state  prediction: 

^<,-1  ~  = 

state  prediction  error  covariance: 


p%k- 1  -  E{[xk  -  *;*_!][**  -  *5fc_ir  w,  z*-1}  =  Fr-i^ri|fc-ii^i + ^-iQr-iGr-i.  (22) 

The  mode-conditioned  predicted  measurement  for  sensor  l  is 

4m’i  :=  br\i^k (23) 

Using  the  linearized  version  (3),  the  covariance  of  the  mode-conditioned  residual 

."MM J(  0  $m,l 

uk  —  zk  ~  zk  > 


is  given  by  (assume  q= 2,  the  case  of  2  sensors) 

sT'1  :=  EWw»rlw'i«r,  z1-') = + *?', 
st’2  :=  £{^'2(i)>'r2<‘,'w,z‘-1) = H73PEk- i«r2' + «r2 
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where  H™’1  is  the  first  order  derivative  (Jacobian  matrix)  of  hm'\.)  evaluated  at  the  state  prediction 
ajfc|fc-i  (see  (23)).  Note  that  (24)  and  (25)  assume  that  z’lf1  originates  from  the  target.  The  results 
(24)  and  (25)  do  not  depend  upon  the  actual  measurements. 

As  mentioned  earlier,  since  our  approach  to  the  problem  deals  not  only  with  the  asynchronous 
measurements  but  also  with  multiple  simultaneous  measurements  [10, 11]  arising  from  two  separate 
sensors  that  are  tracking  a  single  target  through  a  common  surveillance  region,  a  method  for  fusion 
of  multiple  measurements  has  to  be  devised.  In  order  to  do  this,  now  the  combined  covariance  S™ 
of  the  mode-conditioned  residual  obtained  from  (24)  and  (25)  also  needs  to  be  considered  as  follows 


S?  := 


pm  £Vm,l/ 

n b|fc-i  nk 


0  iff 


Step  3.  Measurement  validation  for  sensors  1  and  2  (im  €  Mn)  : 

There  is  uncertainty  regarding  the  measurements’  origins.  Therefore,  we  perform  validation  for 
each  target  separately.  One  sets  up  a  validation  gate  for  sensor  l  centered  at  the  mode-conditioned 
predicted  measurement,  z™'1.  Let  (|A|  =  det(A)) 


ia  :=  arg  {  max  S71’1 1 1 
(m€Mn  K  I J 


Then  measurement  zjj*^  (i=l,2,...,m;)  is  validated  if  and  only  if 

fJ(‘)  ima,l1/rr,m0,!i-1r  l{i)  -maM  ^ 


Zu  -  Z, 


where  7  is  an  appropriate  threshold.  The  volume  of  the  validation  region  with  the  threshold  7  is 
Vl:=cntl  7nW2|5£W|1/2  (28) 

where  nzi  is  the  dimension  of  the  measurement  and  cnzl  is  the  volume  of  the  unit  hypersphere  of 
this  dimension  (ci  =  2,  C2  =  n,  C3  =  47r/3,  etc.).  Choice  of  7  is  discussed  in  more  detail  in  [4,  Sec. 
2.3.2].  After  performing  the  validation  for  each  target  separately,  we  deal  with  all  the  validated 
data  for  measurement  fusion. 

Step  4 •  State  estimation  with  validated  measurement  from  sensors  1  and  2  (im  € 
Mn)  : 

Prom  among  all  the  raw  measurements  from  sensor  l  at  time  k,  i.e.  Zlk  :=  {zkl\zk2\ ..., zkrn^}, 
define  the  set  of  validated  measurement  for  sensor  l  at  time  k  as 

(29) 

where  rhi  is  total  number  of  validated  measurement  for  sensor  l  at  time  k  and 

y2‘):=4(',)  (30) 

where  1<  ii  <  I2  <  ...  <  lmt  <  mi  when  ihi  ^0.  Define  the  association  events  (hypotheses)  as 
follows  (here  we  follow  [10]) 
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•  6^’°:  none  of  the  measurements  in  or  Ffc2  is  target  originated. 

•  only  in  Ffc2  is  a  target  measurement,  all  other  measurements  in  Y£  or  Yj?  are  clutter, 
i  =  0,j  =  1 

•  6]f:  only  yj®  in  F*,1  is  a  target  measurement,  all  other  measurements  in  Y^  or  Y£  are  clutter, 
*  =  l,...,mi,  j  =  0. 

•  @%k  •  yl(l)  yl^  in  Y^  and  Ffc2,  respectively,  are  target  measurements,  all  other  measure¬ 
ments  are  clutter,  i  —  1,  j  =  1, ...,  m2. 

Therefore,  there  are  a  total  of  mim2+mi-t-m,2+l  possible  association  hypotheses,  each  of  which 
has  an  association  probability.  Define  the  mode-conditioned  association  event  probabilities  as 

/3-.U  ;=  P{e)?\M?,Y^YlZk-'}.  (31) 

Exploiting  the  diffuse  model  for  clutter  in  [1,  4],  it  turns  out  that 


where  Pol  and  Pd2  axe  the  detection  probabilities  that  the  sensors  1  and  2  detect  the  target, 
respectively,  Pqx  and  Pc1  are  probabilities  the  target  is  in  the  validation  region  observed  from 
sensors  1  and  2,  respectively,  C  is  a  normalization  constant  such  that  X)S=o  Ph™  =  1  Vm  and 

N[rr;y,P]  :=  |27rP|-1/2exp  {x  -  y)'P_1  (x  -  y)  . 
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Define  the  mode-conditioned  innovations 


as 


9 


where  Kalman  gains,  W ™'1'3 ,  are  computed  as 


=  < 


wm,0,0 

=  0, 

k 

II 

—3 

?r  ^ 

k 

if 

II 

w 

k 

_  pm 

rk\k- 

-  [«r'' 

for  *  =  0,  j  —  0 

yr1'^’1]-1  0],  for  i?0,j  =  0 
>  Hr2'^-2]-1],  for  i  =  0,3^0 
T'IS?}-1,  for  *#0,^0, 

~  f  r  ~  i/  «.  ft/i 

and  H ™  =  \Hk’  Hk  ’  .  Therefore,  mode-conditioned  update  of  the  state  estimate 


(38) 


mi  m 2 


*3*  ~E{xt\MZ‘,Zk-1,Yk1,Yi}  = 

i=0  j=0 

and  covariance  of 

_  mi  777-2  .  .  .  .  .  .  . 

Pk\k  ■=  P%k- !  -  E  E 

«=o,(t,i)#(o,o)  j=o 

i=0j=0 


(39) 


(40) 


mi  7772  .  .  .  . 

£  £  P™™W™'l'3vm™ 

i=0j~0 


777i  m2  .  . 

£  £ 

t=0j=0 


Step  5.  Update  of  mode  probabilities  (Vm  €  Mn)  : 

rt~P[MF\zk]  =  brt-K 


(41) 


where  C  is  a  normalization  constant  such  that  £  p.™  =  1. 

m 

Step  6  Combination  of  the  mode-conditioned  estimates  (Vm  €  Mn)  :  The  final  aug¬ 
mented  state  estimate  update  at  time  k  is  given  by 

**l*  =  £ro*Jb|fcM?  (42) 

and  its  covariance  is  given  by 

Pk\k  =  Y.m  {Pk\k  +  [*fc|k  -  4|fc]  [x™\k  -  xk\k]'^T-  (43) 

From  the  final  augmented  state  (see  (42)),  the  state  filtered  vector  and  the  state  smoothing 
vector  £fc_i|fc  can  be  easily  obtained. 

V  Simulation  Example 


The  following  example  of  tracking  a  highly  maneuvering  target  in  clutter  is  considered.  The  target 
starts  at  location  [21689  10840  40]  in  Cartesian  coordinates  in  meters.  The  initial  velocity  (in  m/s ) 
is  [-8.3  -399.9  0]  and  the  target  stays  at  constant  altitude  with  a  constant  speed  of  400  m/s.  Its 
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trajectory  is  a  straight  line  with  constant  velocity  between  0  and  20s,  a  coordinated  turn  (0.15 
rad/s)  with  constant  acceleration  of  60  m/s2  between  20  and  35s,  a  straight  line  with  constant 
velocity  between  35  and  55s,  a  coordinated  turn  (0.1  rad/s)  with  constant  acceleration  of  40  m/s2 
between  55  and  70s,  and  a  straight  line  with  constant  velocity  between  70  and  90s.  The  target 
motion  models  are  patterned  and  modified  after  [1].  In  each  mode  the  target  dynamics  are  modeled 
in  Cartesian  coordinates  as 

%k  =  +  G^k-  l^r  (44) 

l  +  ^,/e-i^r  (45) 

where  the  augmented  state  of  the  target  consists  of  position,  velocity,  acceleration,  and  the  process 
noise  in  each  of  the  three  Cartesian  coordinates  (a:,  y,  and  z)  at  4  and  4-1  as  well  as  the  delay 
time  dk  at  4.  Thus  both  Xk  and  £kd  are  of  dimension  25  (nx  =  25).  Three  maneuver  models  are 
considered  in  the  following  discussion.  The  system  matrices  Fk,k- 1,  1,  Fkdtk-i  and  Gkdtk-i 

are  defined  as 
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Model  1.  Nearly  constant  velocity  model  with  zero  mean  perturbation  in  acceleration 
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where  T  is  the  sampling  period.  The  standard  deviation  of  the  process  noise  of  M1  is  5  m/s2  (as 
in  [1]). 

Model  2.  Wiener  process  acceleration  (nearly  constant  acceleration  motion) 
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The  standard  deviation  of  the  process  noise  of  M2  is  7.5  m/s2  (as  in  [1]). 

Model  3.  Wiener  process  acceleration  (model  with  large  acceleration  increments,  for  the  onset 
and  termination  of  maneuvers),  with  F3  =  j F2,  G3  —  G2,  Fj  =  F2  and  G\  =  G2d.  The  standard 
deviation  of  the  process  noise  of  M 3  is  40  m/s2  (as  in  [1]). 

The  initial  model  probabilities  are  //q  =  0.8,  /j,2  =  0.1  and  /j.3  =  0.1.  The  mode  switching 
probability  matrix  is  given  by  (as  in  [1]) 
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The  Sensors:  Two  sensors  are  used  to  obtain  the  measurements.  Sensor  1  and  Sensor  2  are 
located  at  |£i,yi,zi]=[-4000  4000  0]  m  and  [a?2 , 3/2 ,  ^2] = [5000  0  0]  m,  respectively,  and  the  central 
processor  is  collocated  with  sensor  1  platform  (we  assume  that  there  is  no  time  delay  between  sensor 
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1  and  central  processor  and  there  is  fixed  but  unknown  time  delay  between  sensor  2  and  central 
processor).  The  measurements  from  sensor  l  for  model  m  are  z[  =  hm’l( xk)  +  w™'1  for  l  =  1  and  2, 
reflecting  range  and  azimuth  angle  for  sensor  1  (radar)  and  azimuth  and  elevation  angles  for  sensor 

2  (infrared).  The  range,  azimuth,  and  elevation  angle  transformations,  respectively,  are  given  by 

n  =  {(x  -  Xif  +  {y  -  yi)2  +  (z  -  zt)2}1/2  (57) 

ai  =  tan-1[(j/-2/i)/(a:-xi)]  (58) 

et  =  t&n~\z-zi)l{{x-xif +  (59) 

As  we  see  from  (1),  (2),  (4)  and  (5),  the  measurements  obtained  from  sensors  1  and  2  can  be 
expressed  as 

z\  =  ^([-^O.O.O.Ojxfc)*^  (60) 

zl  ~  k2([0>0>  ^,fc-i»^,fc-uO]xfc)  +  w%.  (61) 

The  measurement  noise  w™’1  for  sensor  l  is  assumed  to  be  zero-mean  white  Gaussian  with  known 
covariances,  R1  —  diagjqv,  qa{\  =  diag[400m2, 49mrad2]  with  qr  and  qa\  denoting  the  variances 
for  the  radar  range  and  azimuth  measurement  noises,  respectively,  and  R2  =  diag[g02,ge]  = 
diag^mrad2, 47nrad2]  with  qa 2  and  qe  denoting  the  variances  for  the  infrared  sensor  azimuth  and 
elevation  measurement  noises,  respectively.  The  sampling  interval  was  T=ls  and  it  was  assumed 
that  the  probability  of  detection  Pd= 1  for  both  sensors. 

The  Clutter:  For  generating  false  measurements  in  simulations,  the  clutter  was  assumed 
to  be  Poisson  distributed  with  expected  number  of  Aj  =  13  x  10_6/m  mrad  for  sensor  1  and 
A2  =  7  x  10_4/m  mrad  for  sensor  2  [1,  case  1].  These  statistics  were  used  for  generating  the 
clutter  in  all  simulations.  However,  a  nonparametric  clutter  model  was  used  for  implementing  all 
the  algorithms  for  target  tracking. 

Other  Parameters:  The  gates  for  setting  up  the  validation  regions  for  both  the  sensors  were 
based  on  the  threshold  7=  16.  With  the  measurement  vector  of  dimension  2,  this  leads  to  a  gate 
probability  Pq= 0.997  (see  [4,  pages  95-96]). 

Simulation  Results:  The  results  were  obtained  from  100  Monte  Carlo  runs.  Fig.  1  shows  the 
true  trajectory  of  the  target.  Fig.  2  shows  the  delay  estimates  (given  unknown  but  fixed  timing 
mismatch  between  the  two  sensors)  based  on  100  Monte  Carlo  runs.  Fig.  3  shows  the  RMSE  (root 
mean-square  error)  for  the  filtered  state  and  the  smoothed  state  (lag  =1)  in  position,  velocity  and 
acceleration.  It  is  seen  from  Fig.  3  that  the  smoothing  method  shows  better  accuracy  than  the 
filtering  method  as  well  described  in  [9] .  Fig.  4  shows  a  comparison  among  the  performances  of  the 
proposed  IMMMSPDAF  algorithm  dealing  with  asynchronous  measurements  with  unknown  but 
fixed  dk,  with  known  d,  and  the  standard  IMMMSPDAF  algorithm  with  the  assumption  that  d= 0 
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Figure  1:  Trajectory  of  maneuvering  target  (read  left  to  right,  top  to  bottom),  (a)  Position  in  xy 
plane,  (b)  x  and  y  velocities,  (c)  x  and  y  accelerations,  (d)  magnitude  of  accelerations 


always  applies.  It  is  seen  from  Fig.  4  that  when  the  unknown  but  fixed  timing  mismatch  dk  is  more 
than  one  fifth  of  the  sampling  time,  the  performance  improvement  is  significant  compared  with  the 
standard  IMMMSPDAF  algorithm  that  ignores  the  time-delay  d. 

VI  Conclusions 

We  investigated  an  IMMMSPDAF  algorithm  with  asynchronous  measurement  (there  is  unknown 
but  fixed  timing  mismatch  between  sensor  platforms)  for  tracking  a  highly  maneuvering  target  in 
clutter.  The  proposed  algorithm  was  illustrated  via  a  simulation  example  where  it  outperformed 
a  standard  IMMMSPDAF  algorithm  that  ignored  the  possible  timing  mismatch  (especially  when 
the  possible  timing  mismatch  is  more  than  one  fifth  of  the  sampling  time). 
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Figure  2:  Estimation  of  delay  (given  unknown  but  fixed  timing  mismatch  between  two  separated 
sensors)  based  on  100  Monte  Carlo  runs  (read  left  to  right,  top  to  bottom),  (a)  d  =  0.  (b)  d  = 
0.1T.  (c)  d  =  0.3T.  (d)  d  =  0.5T.  (e)  d  =  0.7T.  (f)  d  =  0.9T.  (T  =  sampling  rate) 
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Figure  3:  Comparison  of  filtering  and  smoothing  (lag=l)  for  various  delay  values  (acceleration, 
velocity,  and  position  RMS  errors  (3  rows  each),  read  left  to  right,  top  to  bottom),  (a)  d  =  0.  (b) 
d  =  0.1T.  (c)  d  =  0.3T.  (d)  d  =  0.5T.  (e)  d  =  0.7T.  (f)  d  =  0.9T.  (T  =  sampling  rate).  In  the 
figure  legends,  estimation  refers  to  filtering  and  smoothing  is  with  lag=l. 
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Figure  4:  RMSE  in  position  using  IMMMSPDAF  under  various  scenarios  of  known  delay,  estimated 
delay  and  ignoring  delay,  for  various  delay  values  (read  left  to  right,  top  to  bottom),  (a)  d  =  0. 
(b)  d  =  0.1T.  (c)  d  =  0.3T.  (d)  d  =  0.5T.  (e)  d  =  0.7T.  (f)  d  =  0.9T.  (T  =  sampling  rate).  Unless 
otherwise  stated,  the  results  are  for  filtering. 
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